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Abstract The phase response curve (PRC) is a powerful tool to study the effect of 
a perturbation on the phase of an oscillator, assuming that all the dynamics can be 
explained by the phase variable. However, factors like the rate of convergence to the 
oscillator, strong forcing or high stimulation frequency may invalidate the above as- 
sumption and raise the question of how is the phase variation away from an attractor. 
The concept of isochrons turns out to be crucial to answer this question; from it, we 
have built up Phase Response Functions (PRF) and, in the present paper, we com- 
plete the extension of advancement functions to the transient states by defining the 
Amplitude Response Function (ARF) to control changes in the transversal variables. 
Based on the knowledge of both the PRF and the ARF, we study the case of a pulse- 
train stimulus, and compare the predictions given by the PRC-approach (a ID map) 
to those given by the PRF-ARF-approach (a 2D map); we observe differences up to 
two orders of magnitude in favor of the 2D predictions, especially when the stimu- 
lation frequency is high or the strength of the stimulus is large. We also explore the 
role of hyperbolicity of the limit cycle as well as geometric aspects of the isochrons. 
Summing up, we aim at enlightening the contribution of transient effects in predicting 
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the phase response and showing the limits of the phase reduction approach to prevent 
from falling into wrong predictions in synchronization problems. 

List of Abbreviations 

PRC phase response curve, phase resetting curve. 
PRF phase response function. 
ARF amplitude response function. 

1 Introduction 

The phase response (or resetting) curve (PRC) is frequently used in neuroscience to 
study the effect of a perturbation on the phase of a neuron with oscillatory dynam- 
ics (see surveys in [1-3]). For it to be applied, several conditions are required (weak 
perturbations, long time between stimuli, strong convergence to the limit cycle, etc.) 
so that the system relaxes back to the limit cycle before the next perturbation/kick is 
received. In this case, one can reduce the study to the phase dynamics on the oscilla- 
tory solution (namely, a limit cycle). However, in realistic situations, we may not be 
able to determine whether the system is on an attractor (limit cycle); moreover, the 
system may not show regular spiking, especially because of noise; see for instance [4, 
5]. In addition, even in the absence of noise, strong forcing may send the dynamics 
away from the asymptotic state, eventually close to other nearby invariant manifolds 
[6]; thus, both the rate of convergence to the attractor and the stimulation frequency 
(which can be relatively high; take for instance the case of bursting-like stimuli) play 
an important role in controlling the time spent in the transient state (away from the 
limit cycle). All these factors may prevent the trajectories from relaxing back to the 
limit cycle before the next stimulus arrives and raise the question of the nature of the 
phase variation away from an attractor (that is, in transient states) and how much can 
we rely on the phase reduction (PRC). 

Recently, tools to study the phase variation away from a limit cycle attractor have 
been developed. They rely on the concept of isochrons (manifolds transversal to the 
limit cycle and invariant under time maps given by the flow), introduced by Winfree 
(see [7]) in biological problems, from which one can extend the definition of phase in 
a neighborhood of the limit cycle. In a previous paper [8, 9], we showed how to com- 
pute a parameterization of the isochrons (see also [10-12] for other computational 
methods) as well as the change in phase due to the kicks received when the system is 
approaching the limit cycle but not yet on it. This approach allowed us to control the 
phase advancement away from the limit cycle (that is, in the transient states) and build 
up the Phase Response Functions (PRF), a generalization of the PRCs. In [8], exam- 
ples of neuron oscillators were shown in which the phase advancement was clearly 
different for states sharing the same phase. A review of these tools is presented in 
Sect. 2. 

In Sect. 3, we complete the extension of advancement functions to the transient 
states by defining the Amplitude Response Function (ARF), and we provide methods 
to compute it by controlling the changes induced by perturbations in a transversal 
variable, which represents some distance to the limit cycle. One of the methods pre- 
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sented here to compute the ARFs is an extension of the well-known adjoint method 
for the computation of PRCs; see, for instance, [13, 14] or Chap. 10 in [1]. 

Indeed, the knowledge of both the PRF and the ARF allows us to consider special 
problems in which these functions can forecast the asymptotic phase of an oscillator 
under pulsatile repetitive stimuli. In the case of a pulse-train stimulus, the variations 
of the extended phase and the amplitude can be controlled by means of a 2D map; 
this 2D map extends the classical ID map used when the dynamics is restricted to the 
limit cycle or phase-reduction is assumed; see, for instance Chap. 10 in [1]. Another 
successful strategy to deal with kicks that send the dynamics away from the limit 
cycle is the so-called second-order PRC (see [15-17]), which measures the effects of 
the kick on the next cycle period, taking into account that synaptic input can span two 
cycles. 

As an illustration of the method, in Sect. 4, we then consider a canonical model 
for which we compute the PRFs and ARFs thanks to the exact knowledge of the 
isochrons. In this "canonical" example, we apply a two parametric periodic forcing 
(varying the stimulus strength and frequency) and make predictions both with our 2D 
map and the classical ID map; we use rotation numbers to illustrate the differences 
between the two predictions and we observe differences up to two orders of magni- 
tude in favor of the 2D predictions, especially when the stimulation frequency is high 
or the strength of the stimulus is large. We also use this example to shed light on the 
role of hyperbolicity of the limit cycle as well as geometric aspects of the isochrons 
(see also [18] for a related study of the effect of isochrons' shear). Finally, we also 
present the comparison of the two approaches in a conductance-based neuron model, 
where we do not know the isochrons analytically. 

Summing up, we aim at enlightening the contribution of transient effects in pre- 
dicting the phase response, focusing on the importance of the "degree" of hyperbol- 
icity of the limit cycle, but also on the relative positions of the isochrons with respect 
to the limit cycle. Since PRCs are used for predicting synchronization properties, see 
[19], Chap. 10 in [1] or Chap. 8 in [2], our final goal is to show the limits of the phase 
reduction approach to prevent wrong predictions in synchronization problems. 



2 Set-up of the Problem: Isochrons and Phase Response Functions (PRF) 

In this section, we set up the problem and we review some of the results in [8] that 
serve as a starting point of the study that we present in this paper. 
Consider an autonomous system of ODEs in the plane 

x = X(x), xe£/cR 2 , (1) 

and denote by <p t the flow associated to (1). Assume that X is an analytic vector field 
and that (1) has a hyperbolic limit cycle r of period T, parameterized by 0 = t/T as 

y : T -> R 2 

(2) 

0\-> y(0), 

where T = R/Z, so that y(0) = y(0 + l). 
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Under these conditions, by the stable manifold theorem (see [20]), there exists a 
unique scalar function defined in a neighborhood £2 of the limit cycle r, 

Q : Q C M 2 -> T 

(3) 

xh^0(x) 

such that 

lim \4>t(x)-y(t/T + 0(x))\=O, (4) 

if the limit cycle is attracting. If the limit cycle is repelling, the same is true with 
t -> —oo. 

The value G (x) is the asymptotic phase of x and the isochrons are defined as the 
sets with constant asymptotic phase, that is, the level sets of the function 0. The 
same construction can be extended to limit cycles in higher dimensional spaces, but 
since the applications in this paper will be restricted to planar systems, we give the 
definitions in R 2 . 

Moreover, we know from [21] that there exists an analytic local diffeomorphism 
^:Tx [cr-,cr+] M 2 

(5) 

(<9,<r)i-> K(0,a), 

satisfying the invariance equation 

where T is the period and X is the characteristic exponent of the periodic orbit. 

We can describe (6) as saying that if we perform the change of variables given by 
K, the dynamics of the system (1) in the coordinates (0, a) consist of a rigid rotation 
with constant velocity 1/ T fox 6 and a contraction (if A < 0) with exponential rate 
X/T for a. That is, 

0 = 1/T, 

( 7 ) 

a = Xcr/T, 

and </>t(K(0, a)) = K(0 +t/T, cre (k/T)t ), where fa is the flow associated to (1). See 
Fig. 1 for an illustration of the evolution of the variables (0, a) along an orbit of the 
system. 

Remark 2.1 In [8], we presented a computational method to compute the parameter- 
ization K defined in (6) numerically. 

Given an analytic local diffeomorphism K satisfying (6), we know from Theo- 
rem 3.1 in [8] that the isochrons are the orbits of a vector field Y satisfying the Lie 
symmetry [Y, X] = fiY with /jl = X/T. That is, 

Y oK(0,cr) = d a K(0,cr). (8) 
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Fig. 1 The new coordinate 
system. Representation of a 
limit cycle (red) and an isochron 
(blue) corresponding to the 
phase level sets 9 = 6q. The 
black dashed curve is a piece of 
the trajectory of the vector field 
X starting at a generic point 
K(0 0 ,a 0 ) 



K(9 0 ,a 0 ) 



K(9 0 ,croe x 



r=W)} = {K(e 1 o)} 



Let us assume that a pulse of small modulus € instantaneously displaces the tra- 
jectory in a direction given by the unit vector w. Mathematically, we consider 

x = X(x) + €w8(t-t s ), (9) 

where € and 8(t) is the Dirac delta function. Then we define the phase response 
function (PRF) as the infinitesimal rate of change of the phase in the direction w of 
the perturbation, that is, 

PRF(x) = D w 0(x) = (w, V<9(x)), 

where Q is the phase function defined in (3) and (•, •) denotes the dot product. In [8], 
we showed that 



V&(K(0,cr)) = . 

V n T((d a K)^,d e K) 

where (3 a K) 1 - = J (d G K) and the matrix / is given by 



(10) 



'-(! -;)■ 

We will use this notation for the rest of the paper. 

In neuron models, the usual situation is x = (V, . . . ) and w = (1, 0, . . . , 0); thus, 
the phase response function (PRF) is defined as 

PRF(x) = 3y(9(x), (12) 

where 3y denotes partial derivative with respect to the variable V. 



3 The Amplitude Response Function (ARF) 

A pulse stimulation displaces the trajectory away from the limit cycle, producing a 
change both in the phase 0 and the transversal variable a, that we will refer to as 
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the amplitude variable. In our notation, the amplitude variable is a distance measure 
defined by the time from the limit cycle along the orbits of the auxiliary vector field 
Y, transversal to X, defined in (8). In fact, the orbits of the vector field Y are the 
isochrons; see, for instance, the blue curve in Fig. 1. The phase-reduction approach 
assumes that the amplitude decreases to zero before the next pulse arrives and, there- 
fore, the amplitude is always zero at the stimulation time. However, if one wants to 
consider pulses that arrive before the trajectory relaxes back to the limit cycle, one 
needs to compute also the amplitude displacement in order to predict the coordinates 
of the point at the next stimulation time. 

In this section, we introduce the amplitude function and Amplitude Response 
Function (ARF), the analogues of the phase function (3) and the PRF (12) for the 
variable o . Finally, we provide a formula to compute them given the diffeomorphism 
K introduced in (5). 

3.1 Definitions 

Given an analytic local diffeomorphism K as in (5) satisfying (6), it follows that there 
exists a unique function U, defined in a neighborhood Q of the limit cycle r, 

H : Q c M 2 -> R 

(13) 

xi-> 27 (x) 

such that 

Z((t) t (x)) = Z(x)e kt / T , 

where <p t is the flow associated to the vector field X. The level curves of U are closed 
curves that we will call amplitude level curves or, in short, A-curves. 

Analogous to the phase isochrons, it can be seen that given an analytic local dif- 
feomorphism K, as in (5), satisfying (6), the A-curves are the orbits of a vector field 
Z, satisfying [X, Z] = [Z,X]= 0; see the Appendix for a proof of this result. More 
specifically, 

Z(K(P,cr)) = d 9 K(P 9 cr). (14) 

Expressed in the variables (0, a) introduced in (5), the motion generated by Z is 
given by {0 = 1, a = 0}. 

A pulsatile kick in the direction given by the unit vector w (see (9)) causes a 
change in the amplitude variable. Analogous to the PRF introduced in (10), we define 
the Amplitude Response Function (ARF) as the infinitesimal rate of change of the 
phase in the direction w of the perturbation, that is, 

ARF(x) = D w E(x) = (w, V27(x)). 

In neuron models, the ARF typically measures the change in amplitude under the 
action of a pulsatile kick in the direction of the voltage V, that is, 

ARF(x) = a y £(x), 

where 3y denotes partial derivative with respect to the variable V. 
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3.2 Computation of the PRFs and the ARFs 

In this section, we provide a formula to compute the functions V<9 and VZ* given 
the diffeomorphism K introduced in (5). 

Using the parameterization K introduced in (5) and writing K(x, y) = K(0 (x, y), 
U(x, y)) = (K x , K y ), where 0 and U are the functions introduced in (3) and (13), 
respectively, we have that 



5i K x 9cr K x 



d x 0 dy@\ _ (\ 0 



o 1 



and, therefore, V<9 = (d x 0, d y 0) and VI7 = (d x E, d y U) are given by 

-l 



3# ^ x 3 a 

36) A'-y 3 a 

1 (dgKy 

(d a K^,d e K) \-d e K y 



-d*K x 



Hence, 



V6)(^(6>,cr)) = 
V 17(^(0, or)) = 



3^(0,0-) 



(d a KH0,cr),deK(G,o)Y 
(d a K^(0,a),d e K(0,cr))' 



and 



(15) 



Using the vector field description given in (8) and (14), we can rewrite the expres- 
sion above, using the relations d a K = Y o K and d$K = Z o K, as 



V<9(£(0,cr)) = — : 



K(0,cr) 



and VE(K(0,cr)) = 



(Z ± ,Y) 



K(0,cr) 



By the invariance equation (6), we know that X = jZ + jgY and, therefore, 



V0(K(O,a)) = 



T(Y ± ,X) 
Xa Z L 



and 



K(p,cr) 



T (Z L ,X) 



(16) 



K(6,cr) 



Remark 3.1 Notice that expression for VU in (16) might suggest that it has a singu- 
larity at a = 0. Nevertheless, the vanishing terms in the numerator and denominator 
cancel out at a = 0, and using that Z(K(0, 0)) = d e K(0, 0) = X(K(0, 0)), the value 
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at o = 0 is given by 

X ± (y(0)) 



V27(y(0)): 



where tfi(0) = 3 (J ^(6>, 0) = Y(K(0, 0)). 
3.3 The Adjoint Method for the ARF 

The most common method to compute the PRC, the so-called adjoint method, uses 
that the function V<9 evaluated on the limit cycle r is a periodic solution of some 
adjoint equation (see, for instance, [1]). In the generalization introduced in [8], it 
was shown that the adjoint method can be extended to compute V<9 for points in a 
neighborhood of the limit cycle, for which the periodicity condition is not satisfied. 
Indeed, Q = V0 satisfies the equation 

^ = -DX T (Mp))Q, (17) 

where DX T is the transpose of the real matrix DX. In this case, the method just 
requires an initial condition, so that it can be solved uniquely. The initial condition is 
provided by formula (15). 

The same result can be extended to compute the change in the transversal variable 
a due to a pulse stimulation. In the following proposition, we provide the differential 
equation satisfied by VE(p) where p = K(0, a) is a point in a neighborhood Q of 
the limit cycle y evolving under the flow of X. 

Proposition 3.2 Let r be a hyperbolic T -periodic orbit of a planar analytic vector 
field X parameterized by 0 according to (2). Assume that there exists a change of 
coordinates K in a neighborhood Q satisfying (6). Then the function WE along the 
orbits of the vector field X satisfies the adjoint equation 



dQ (X 
dt 



(^-DX T (MP)))Q, (18) 



where <p t is the flow of the vector field X and A is the characteristic multiplier of the 
periodic orbit, with the initial condition 

_ kZ(p) Z ± (p) 

2(0) = t (zHp),x( P ))' (19) 

where Z ± (K(6, a)) = Jd e K(0, a). 

Proof We will show that the function V E evaluated along the orbits <fi t (p) of X 
satisfies the adjoint equation (18). From expression (16), we have that 

1 ' t (zH<t>t(p)),x(4> t (p))) 
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We now compute the derivative of WU(4> t (p)) with respect to time. In order to 
simplify notation, we set x := (p t (p)- We will also use that Z 1 - = JZ where / is the 
matrix (11). Using that f t Z(x) = DZ(x)Z(x), we have from (20) 

d WE(x) _ X(dZ/dt)(x)JZ(x) + XZ(x)JDZ(x)X(x) 
dt ~ T(Z ± (x),X(x)) 

XE(x)JZ(x)({JDZ(x)X(x), X(x)> + (/Z(x), DX(x)X(x))) 
T(ZH*),X(x)) 2 ' 

Using that DXZ = DZX, expression (20) and dot product properties (namely, 
(JZ(x), DX(x)X(x)} = (DX(x) T JZ(x), X(x)», we obtain 

d (X I T) X E (x) J Z (x) + X U (x) / D X (x) Z (x) 

— VI7(x) = 

dt Tg(x) 

VU(x)(JDX(x)Z(x) + DX(x) T JZ(x), X(x)) 
(Z^(x),X(x)) ' 

Using JDXix) + DX(x) T J = trace(DX)(x)/, and denoting r(x) := trace(£>X)(x), 
we are led to 

d_ = (X/T - Z)X(x) r + t(x))XE(x)JZ(x) _ W E(x)((t(x)JZ(x), X(x))) 

dt r(Z- L (x),X(x)) (Z- L (x),X(x)> 

Finally, using again (20), we have 

— V17(x) = (A/r - Z)X(x) r + t(x))V17(x) - V27(x)r(x) 

dt v 7 

= (X/T - DX(x) T )WE(x), 
as we wanted to prove. □ 

4 Periodic Pulse-Train Stimuli 

The purpose of this section is to show the convenience of using the response func- 
tions away from the limit cycle to obtain accurate predictions of the ultimate phase 
advancement. To this end, we force a system with pulse-trains of period 7^ <^ To for 
trajectories near a limit cycle r of period 7b and characteristic exponent X. 

Given an oscillator, assume that it is perturbed with an external instantaneous stim- 
ulus of amplitude € in the voltage direction every T s time units, that is, 

N 



i = Z(x) + €W^«(r-jT J ), (21) 

7=0 



where w=(l,0),€^Cl and 8 is the Dirac delta function. This system can represent, 
for example, a neuron which receives an idealized synaptic input from other neurons. 

Springer 



Page 10 of 26 



O. Castejon et al. 



Remark 4.1 In the sequel, we will also use co s = 1/T S , the frequency of the stimulus, 
and o)o = l/ To, the frequency of the limit cycle F. Then the quotient co s /coo indicates 
how many inputs the oscillator receives in one period. 

In order to know the evolution of the perturbed oscillator after each time period T s , 
it is enough to know how the variables 6 and a change. We recall that the variation of 
the variable 6 produced by an external stimulus is given, to first order in the stimulus 
strength € , by the PRE Similarly, the variation of the variable o is given to first order 
by the ARF. Hence, we can consider the following map, which approximates the 
position of the oscillator at the moment of the next kick: 

0 n + l =0 n +€ PRF(0„, O n ) + J (mod 1), 

T o (22) 
cfn+i = (°n + € ARF(6 n , a n ))e kTs / T °. 

Moreover, we can compare it with the map obtained by considering the classical PRC 
(see, for instance, Chap. 10 in [1]), which is 

e n +l=e n +€¥RC(e n ) + ^- (modi). (23) 
To 

In the latter case, we are assuming that the perturbation happens always on the limit 
cycle and, therefore, a n = 0 for all n. The possibility that this might not be a realistic 
assumption (for example, if the stimulus period T s is too small, the limit cycle is 
weakly hyperbolic or the strength of the stimulus € is too large) has been already 
pointed out in the literature; see, for instance, [22] or Chap. 10 in [1]. However, 
other factors could play a role, for example, the geometry of the isochrons (curvature, 
trans versality to the limit cycle, etc.). Our aim is to consider some examples and see 
in which cases the ID map (23) gives a correct prediction or, on the contrary, one 
requires the 2D map (22) to correctly assess the true dynamics of the phase variable. 

To quantify the long-term agreement or disagreement between the ID and the 2D 
predictions, we compute an approximation of the rotation numbers after N iterations 
for both maps. More precisely, given an initial condition (Go, &o), we compute 

1 N 

P= lim TT T\{9j-0j-i), (24) 

7 = 1 

where 6 denotes the lift of 0 to R. Then, for the 2-dimensional map (22) and assuming 
N large enough, the rotation number can be approximated by 

T 1 

P2D:=^ + -^PRF(^,a ; ), (25) 
To N ^ 



and by 



T 1 N ~ l 

Piz) := ^ + -e y>RC(0;), (26) 
T ° N 7=0 



for the 1 -dimensional map (23). 
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These approximate rotation numbers will be our main indicator to compare the 
dynamics predicted by the ID map with that of the 2D map. In order to dissect the 
causes that create the eventual differences between the two maps and highlight the 
shortcomings of the phase-reduction approach, we have first considered a "canonical" 
example in which the isochrons can be computed analytically. Next, we consider a 
conductance-based model, in which the isochrons have to be computed numerically 
and we obtain similar comparative results. 

4.1 Examples 

4.7.7 A Simple Canonical Model 

We consider a simple model having a limit cycle with two parameters, a and a, 
that control the hyperbolicity of the limit cycle and the angle between the isochrons 
and the limit cycle, respectively. The system has the following expression in polar 
coordinates: 



r = ar{\ — r 2 ), 
cp = 1 -f aar 2 , 

for a, a e R, and the following one in Cartesian coordinates: 

x =ax(l - (x 2 + y 2 )) - y(l +aa(i 2 + y 2 )), 
y = ay(l - (x 2 + y 2 )) + x(l + aa(x 2 + y 2 )). 



(27) 



(28) 



The limit cycle corresponds to r = 1 and the dynamics on it are given by <p = 
1 + aa\ therefore, cp(t) = cpo + (1 + aa)t mod2jr . The period of the limit cycle r 
is 7b = 2tt/(1 + ad). A parameterization of the limit cycle in terms of the phase 
0 = t/T 0 , for 0 e [0, 1) is y(0) = (cos(2tt(9), sin(27T(9)). 

Now, consider the vector field F, given in the polar and Cartesian coordinates by 

r =ar 3 , and i = a (x 2 + y 2 ) (x + ay) , 
cp = —aar 2 ; y = a(x 2 + y 2 )(y — ax). 

It is easy to check that Y satisfies [X, Y] = — 2aY. Then, using (8), we find that 
l± = —2a and 



KUKcr) = j J- — l - — cos^2tt(9 + ^a\n{\-2aa) 



' - <- • 1 n. 



sin 2tt6 + -a\n(l -2ao) , (29) 



1 — 2ao \ 2 J J 

with 6 £ [0, 1) and a > l/(2a). 
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Notice that the function K can be easily inverted using that r 2 = x 2 + y 2 = 
(1 - 2aa)~ l and arctan(^) = 2it0 + \a ln(l - 2aa). Then K~ l (x, y) = (<9(jc, y), 
£(x,y)), where 



«9(x,y) = i-(arctang)-I fl ln(l)), 

Thus, the dynamics for (0, a) are given by 

# = 1/ 7b, 
<T = — 2aa. 



1 / 1 



(30) 



The vector field Z, defined in (14), has the following expression in Cartesian co- 
ordinates and polar coordinates, respectively: 

x — —2Tty, and r = 0, 
y = 2nx; (p = 2it. 

Therefore, we find that VG(p) = ^2 ( — ^ +ax,x + ay), and, by the parameter- 
ization y of the limit cycle, V<9(y(0)) = ^(— sin(2jr^) + <2 cos(27r#), cos(27r#) + 
asin(27i6>)). Similarly, VU(p) = (^4,^4), and VZ(y(0)) = (cos(2tt(9), 
sin(2jr^)). 

From the last equations, we can then obtain 

PRF(^(6>,a)) = -"^-^ ( sinf 2tt(9 + -aln(l - 2aa) J 

— acos\2jr0 + -<zln(l — 2aa) J J (31) 



and 



(I — 2aa) 5/2 ( 1 \ 

ARF(^(6>,a)) = — - — cos( 2tc6 + -a ln(l - 2aa) j . 

In Fig. 2, we show the PRF and the ARF for a specific isochron for representative 
values of the parameters, a = 10 and a = 0.1. An important remark is that the PRF 
is far from being constant along isochrons, whereas the ARF is clearly nonzero near 
the limit cycle (cr = 0). These features will have a significant effect when comparing 
the predictions provided by the ID map and the 2D map. 

We want to stress the role of the parameters a and a. On one hand, the parameter 
a determines the hyperbolicity of the limit cycle, since its characteristic exponent is 

A, = —2a To. 

Hence, for small values of a the contraction to the limit cycle will be weak, while as 
a goes to infinity A tends to — Ait I a. On the other hand, the parameter a determines 
the transversality of the isochrons to the limit cycle. Indeed, denoting by the angle 
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(a) a = 0.1 (b)a = 10 

Fig. 3 Limit cycle and isochrons. The limit cycle (red) of system (28) and some isochrons (blue) for 
different values of the parameter a. In both cases, a = 10 



between the isochron {p eM 2 : 0{p) = 0} and the limit cycle at the point y(0), we 
have 

(y'(Q),v®Hy(Q))) 

cos 6 = : . 

Wy'(O)\\\\v0Hy(P))\\ 

Computing explicitly the right-hand side of the equality, it is straightforward to verify 
that 

lira 

cos 6 = 



In particular, note that ft is independent of the variable 0 . Moreover, for a = 0 the 
isochrons are orthogonal to the limit cycle and, as a tends to infinity, they become to 
it (see Fig. 3). 
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4.7.2 Numerical Simulations 

In this section, we use the analytic expressions obtained in the previous subsection 
for the PRF, ARF, and PRC to compute and compare the maps defined in (22) and 
(23). Moreover, as we also have an explicit formula for the parameterization K and 
its inverse ^ _1 ,we can integrate numerically system (28), perturb it periodically, and 
obtain a sequence (x n , y n ). Then we can compute analytically 

(O n ,a n ) = K- 1 (x n ,y n ), (32) 

and compare it with the iterations obtained using the maps (22) and (23). In the 
following, we will call the approximation of the rotation number obtained by this 
method simply p, to distinguish it from p2D and p\o defined previously in (25) and 
(26), respectively. The following lemma gives a description of the dynamics expected 
in the 1 -dimensional map. 

Lemma 4.2 For k e Z, let us denote 

c ' = -(f- 

Then, the fixed points of the 1-dimensional map (23) can be computed analytically 
and: 

• If 1 + a 2 + C\ < Ofor all k e Z, the map (23) has no fixed points. 

• If there exists k e Z such that 1 + a 2 + C\ < 0 and 



aC k + Jl+a 2 -C 2 



l+a 2 

the map (23) has the fixed point 



< 1, 



2n 

Moreover, ifa<Ck and 



I / —aCjc + J I + a 2 — C 

0 1 = — — arccos 



l+a 2 

there exists also another fixed point 



aC k -Jl+a 2 -C 2 

< 1, 



I /-aC k -^\+a 2 -C 2 
= — arccos I 



2tt \ \-\-a 1 



Proof The fixed points 6* of map (23) must satisfy 

ePRC(<9*) + — =k 



4y Springer 



Journal of Mathematical Neuroscience (2013) 3:13 



Page 15 of 26 



for some k e Z, or equivalently 

PRC(<9*) + — =0. 
v 7 lit 

Substituting PRC(#*) in the above equation by expression (31) with a = 0 and rear- 
ranging terms we have 

sin(27r(9*) = a cos(2tt(9*) + C k . (33) 

Taking squares of both sides of the equality and using trigonometric properties, we 
obtain 

(1 + a 2 ) cos 2 (2tt(9*) + 2aC k cos(2tt(9*) + c\ - 1 = 0, 

which is an equation of degree 2 in cos(27r#*). Solving it, after some simplifications, 
we obtain 



-aC k ± Jl+a 2 -C 2 
l+a 

It is clear that for Eq. (34) to have real solutions, the right-hand side must have mod- 



cos (2jt(9*) = f-—^ . (34) 



ulus at most 1 and 1 + a 2 — C& > 0. In this case, the solutions of (34) are 



I / —aCk ± J 1 + a 2 — C 



arccos 



In \ l+a 2 



However, as we have taken squares in Eq. (33), we still have to check whether 6+ and 
0* are solutions of (33). It is easy to check that 0* always solves (33), while 6- is a 
solution only when a < . □ 

Remark 4.3 A natural question is whether the 2-dimensional map (22) and the se- 
quence (32) have the same qualitative behavior. As an example, let us take € = 0.03, 
a = 0.1 and a = 10. In this case, there exists just the fixed point 0+ for the 1- 
dimensional map (23). So, let us take the initial condition (Oo,cro) = (#+>0) an d 
compute its iterates by the three different maps (22), (23), and (32). In Fig. 4, we 
plot the sequences K(0 n ,cr n ) (for clarity, we have just plotted those with n > 200). 
As one can see, map (23) fails to predict correctly the qualitative behavior of the so- 
lution, since (32) seems to be attracted to a quasi-periodic orbit and not a fixed point. 
On the contrary, map (22) correctly predicts this qualitative behavior. 

From now on, we will take the initial condition to be (#o, tfo) = (0.8, 0), that is, 
(*0> yo) ~ (0.30901, —0.95106). In order to explore the effect of both the hyperbolic- 
ity and the transversality of the isochrons to the limit cycle, we will plot the different 
approximations of the rotation numbers p, p2D, and pm for different values of the 
parameters a and a . 

First of all, we will take a = 0.1 and a = 10. This corresponds to considering a 
weakly hyperbolic limit cycle with isochrons that are almost tangent to it. In Fig. 5, 
we show the rotation numbers obtained for different amplitudes and for two different 
stimulus periods 7^. In this case, in order to make the rotation number pm stabilize, 
we have taken N = 1000. 
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Observe also the agreement with the result in Lemma 4.2, which predicts the ap- 
pearance of the fixed point of the ID map when 1 + a 2 — C\ — 0, that is, when 
C| = 101 or, equivalently after substituting T s = To/m, e = 27r/(Vl01m). The fixed 
point appears at € ~ 0.0125 for m = 50 (panel (a) in Fig. 5) and € ~ 0.0312 for 
m = 20 (panel (b) in Fig. 5); both values coincide with the downstroke of the corre- 
sponding values of pip. 

One can see that the rotation number obtained with the 1 -dimensional map 
diverges from the analytical computation, while the one obtained with the 2- 
dimensional map does not. This wrong prediction by the 1 -dimensional approach 
is consistent for all intermediate values of 7^ (not shown here). We point out that, 
although the difference between the 1 -dimensional approach and the other two seems 




Fig. 5 Rotation numbers as a function of the stimulus strength. Rotation numbers as a function of the stim- 
ulus strength for parameter values a = 0.1 and a = 10. Stimulation periods are aT s = 0.0628319 ~ Tq/ 50, 
bT s =0.1570800^ T 0 /20 
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Fig. 6 Iterates for a weakly 
hyperbolic limit cycle. First 100 
iterates of sequences K(6 n ,a n ) 
computed using the three 
different methods. Parameter 
values are taken to be a = 0. 1 , 
a = 10, 6 = 0.022 and 
T s =0.0628319 ^T 0 /50 



Limit cycle - 
1D map 
2D map 

++_ + +^ + + + + ++ Analytic computation 




Fig. 7 Differences among the 
2D prediction and the analytic 
rotation numbers. Absolute 
difference between the rotation 
number obtained with the 
2-dimensional approach and the 
analytic one, that is \pid ~ PV 
in the two -parametric space 
Orel><0 




a> re| , stimulus relative frequency 



rather small (it ranges from 10 -3 to 10 -2 ), we can identify a wrong prediction of 
the qualitative behavior of the system by the 1 -dimensional map. Indeed, in the cases 
where p\D ~ 0 but P2D, P^0, the 1 -dimensional map (23) has a fixed point, while 
the other two do not (see Remark 4.3). For example, in Fig. 6, we plot in the phase 
space the first 100 iterates of the sequences K(6 n ,cr n ), where (0 n ,cr n ) are obtained, 
respectively, using the 2-dimensional map (22), the 1 -dimensional map (23), and ex- 
pression (32). While for the ID map a fixed point is reached, for the 2D and the 
analytic approaches it seems that the dynamics are not so simple. Observe that this 
different qualitative behavior is obtained in spite of the initial condition being on the 
limit cycle. 

Another visualization of the agreement or disagreement between the different ap- 
proximations of the rotation numbers is provided in Figs. 7, 8, and 9. We show the 
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Fig. 8 Differences among the 
ID prediction and the analytic 
rotation numbers. Absolute 
difference between the rotation 
number obtained with the 
1 -dimensional approach 
(phase-reduction hypothesis) 
and the analytic one, that is 
\p 1D -p|,inthe 
two-parametric space (co Ye \, e) 




20 23 26 29 32 35 38 41 44 47 50 



Fig. 9 Ratio between errors 




0.04 


given by the 2D prediction and 




0.0364 


by the ID prediction. Ratio of 




the absolute difference between 




0.0328 


the 2-dimensional approach and 


-C 
O) 


0.0292 


the analytic one (numerator, see 


c 

CD 


0.0256 


Fig. 7) over the absolute 


CO 




difference between the 


00 


0.022 


1 -dimensional approach and the 


£ 


0.0184 


analytic one (denominator, see 


CO 


0.0148 


Fig. 8), that is, 






\P2D ~P\/\P1D ~ P\ 




0.0112 






0.0076 






0.004 



26 29 32 35 38 41 44 47 

co , stimulus relative frequency 



I 



differences between them depending on both € (that is, the strength of the stimu- 
lus) and &> re i := cd s /ooo = To/ T s (the ratio between the frequency of the stimulus 
and the frequency of the limit cycle). In Fig. 7, we plot the absolute difference be- 
tween the rotation number obtained with the 2-dimensional approach and the analytic 
one, namely \p2D — p\, whereas in Fig. 8, we plot the error when using the phase- 
reduction hypothesis, namely \p\o — p\. Both errors are compared in Fig. 9, where 
the ratio \p2D — p\/\piD — P I is displayed. As expected, one can see in Figs. 7 and 8 
that, fixing the stimulus period T s , the worst approximations of p given respectively 
by P2D and p\u are obtained for high values of e. However, fixing the strength of the 
stimulus €, the results for both cases are different: while for the 2-dimensional map 
the worst results are for high frequency ratios &> re i, for the ID approach the worst 
results are obtained, in general, for low co rQ \. Finally, as we also expected, in Fig. 9 
we can appreciate that the 2D approach is always better than the ID. Moreover, the 
difference between p2D and p is, in the worst case, two orders of magnitude smaller 
than the difference between p\D and p. 

As we mentioned above, we use this example to help us understanding the effect 
of the hyperbolicity of the limit cycle and the transversality of the isochrons to it 
in the validity of the PRC approach. In Figs. 10, 11, and 12, we plot the different 
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Fig. 10 Effect of weak hyperbolicity and normal isochrons in the phase prediction. Rotation numbers 
for different stimulus strengths in case of weak hyperbolicity and normal isochrons (a = 0.1 and a = 0). 
Stimulation periods are aT s = 0.0628319 % T 0 /50, b T s = 0.1570800 « 7b/20 




Fig. 11 Effect of strong hyperbolicity and almost tangent isochrons in the phase prediction. Rotation 
numbers for different stimulus strengths in case of strong hyperbolicity and almost tangent isochrons 
(a = 10 and a = 10). Stimulation periods are aT s = 0.0628319 % 7q/50, b T s = 0.1570800 % r 0 /20 



approximations of the rotation numbers varying the parameters a (a = 0 meaning 
loss of hyperbolicity) and a (a = 0 meaning isochrons normal to the limit cycle). 
On one hand, when the limit cycle is strongly hyperbolic (for instance, a = 10 as in 
Figs. 11 and 12), all approximations give a very similar result. Hence, in these two 
cases (even when the isochrons are almost tangent to the limit cycle, which corre- 
sponds to Fig. 11), the use of PRFs and ARFs instead of PRCs seems not necessary. 
In fact, that is what one can expect intuitively: if the attraction to the limit cycle is 
very strong, the system relaxes back to the asymptotic state very quickly, so that at 
each kick we can assume that the state variables are on the limit cycle. Of course, this 
will depend also on the frequency of stimulation co s . 

On the other hand, in Fig. 10, where the contraction to the limit cycle is slow but 
the isochrons are almost orthogonal to the limit cycle, one can see that the ID ap- 
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Fig. 12 Effect of strong hyperbolicity and normal isochrons in the phase prediction. Rotation numbers 
for different stimulus strengths in case of strong hyperbolicity and normal isochrons (a = 10 and a = 0). 
Stimulation periods are a T s = 0.0628319 % T 0 /50, b T s = 0.1570800 % T 0 /20 



proach diverges from the 2D approach and the analytic one. However, for the range 
of € and the two different stimulation periods T s (panels (a) and (b)) considered in 
Fig. 10, the ID prediction still gives a fairly good approximation. Moreover, unlike 
the case where a = 0.1 and a = 10 (see Fig. 5), the ID approach predicts a sim- 
ilar qualitative behavior as the other two approaches. The results for € = 0.04 in 
Fig. 10(a) raise another interesting question since the analytic rotation number p sud- 
denly diverges from the ID and the 2D rotation numbers. This is due to the fact that 
the iterates of the analytic map suddenly fail to encircle the critical point of the con- 
tinuous system (located inside the limit cycle) while the iterates of the ID and the 2D 
maps still do it. Thus, the rotation number for the analytic case may not give accurate 
information. 

In conclusion, it seems that for the 2D map to represent a qualitative improvement 
with respect to the ID it is necessary to have the combination of weak hyperbolicity 
of the limit cycle and "weak trans versality" of isochrons to it. However, the role 
of hyperbolicity seems to be much more important, since in the presence of strong 
hyperbolicity the use of the 2D approach seems completely unnecessary, but for weak 
hyperbolicity the differences between the ID and the 2D maps are present also when 
the isochrons are orthogonal to the limit cycle. 

Remark 4.4 Of course, considering a stimulus strength € large enough, both maps 
(22) and (23) will not give correct predictions, since they are based on first-order 
approximations. In this case, one should consider PRFs of second (or higher) order 
to obtain a correct result; see, for instance, [23, 24] for higher-order PRCs. One has 
to distinguish between these higher order response functions in terms of the stimulus 
strength from the second-order PRCs above mentioned (see [15] for instance) that 
relate to the second cycle after the stimulus. 

In the next example, we apply the same methodology to a more biologically in- 
spired case: a conductance-based model for a point-neuron with two types of ionic 
channels. 
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Fig. 13 The limit cycle of the 
conductance-based model and 
its isochrons. The limit cycle 
(red) and equally spaced in 
phase isochrons (blue) for 
system (35) and 7 app = 190 




4.1.3 A Conductance-Based Model 

We consider a reduced Hodgkin-Huxley-like system, with sodium and potassium 
currents, and only one gating variable: 

V = —^(gmmoo(V)(V - V m ) + g K n(V - V K )+g L (V - V L ) - / app ), 

C m (35) 

h = n 00 (V) -n, 

where V represents the membrane potential, in mV, n is a nondimensional gating 
variable and the open-state probability functions are 



m 00 (V) = 
noo(V) = 



1 



1 + exp(-(V - Vmax,m)/£m) 
1 

1 + exp(-(V - V m3X , n )/k n ) ' 



The parameters of the system are C m = 1 |xF/cm 2 , gNa = 20mS/cm 2 , VNa = 
60 mV, g K = 10 mS/cm 2 , V K = -90 mV, g L = 8 mS/cm 2 , v L = -80 mV, V max , m = 
—20 mV, k m = 15, Vmax,n = —25 mV, k n = 5. 

Here, we will take 4 pp = 190 |xA/cm 2 . In this case, the system has a limit cycle 
with period 7b ~ 1.3055442, and its characteristic exponent is X ~ —0.6055956. That 
is, the limit cycle is weakly hyperbolic, and hence we expect that the 2-dimensional 
approach will give qualitatively different results with respect to the 1 -dimensional 
approach. Figure 13 shows the limit cycle and its isochrons. 

In Fig. 14, we show the PRF and the ARF on the limit cycle (cr =0), panels (a) 
and (b), and for a specific isochron (0=0), panels (c) and (d). 

Remark 4.5 We have chosen a value of the parameter 7 app = 190 for which the system 
presents weak attraction to the limit cycle. However, for this value of 7 app , system (35) 
is not a model of a spiking neuron, but one with high voltage oscillations. Thus, this 
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Fig. 14 Phase and amplitude response functions along isochrons and a -level sets. Phase response function 
(PRF) and amplitude response function (ARF) for system (35) and / app = 190 on the limit cycle (cr = 0), 
panels a and b, and for a specific isochron (6 = 0), panels c and d 



example is not intended to deal with a realistic setting of spike synchronization, but 
to illustrate how to deal with the tools introduced in this paper in the case where one 
does not explicitly have the parameterization K. 

Remark 4.6 In order to compute the parameterization K and the PRFs we have used 
the methods proposed in [8]. The same ideas can be applied to compute the ARFs. 
Briefly, the method consists of two steps. First, to compute the value of a given ARF 
near the limit cycle, where the numeric approximation of the parameterization K is 
valid, expression (16) is used. Second, to compute the value of some ARF far from 
the limit cycle, we just integrate the adjoint system (18) backwards in time using an 
initial condition for ARF close to the limit cycle. 

Again, we have computed the rotation numbers as defined in (25) and (26) varying 
the strength of the stimulus € with fixed stimulation periods. We have taken N = 100 
and initial conditions Go = 0.089 and oo = 0. The results, for two different stimuli 
periods T s , are shown in Fig. 15. Again, note that although the dynamics begin on 
the limit cycle (since ao = 0), the behavior of the 1 -dimensional approach and the 2- 
dimensional approach are quite different. Moreover, for e > 0.4 we find that p\o ~ 0, 
while p2D ~ 0.02. This can be interpreted, similarly to the previous example, as an 
indicator that the ID map (23) has a fixed point, while the 2D map does not. Further- 
more, this indicates that after 100 iterations of the 2D map (22), the state variables 
have turned approximately twice around the fixed point, as one can see from the plots 
of the sequences K(0 n , a n ) computed using both maps (see Fig. 16). 
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Fig. 15 Rotation numbers for different stimulus strengths in the conductance-based model. Rotation num- 
bers for different stimulus strengths and fixed stimulus periods for system (35) 



Fig. 16 Iterates for a weakly 
hyperbolic limit cycle of the 
conductance-based model. 
Sequences K(6 n , a n ) computed 
using the 2-dimensional map 

(22) and the 1 -dimensional map 

(23) , respectively, for system 
(35). The strength of the 
stimulus is e = 0.574604, while 
the stimulation period is 

T s =0.026111 ^r 0 /50 




-24 -22 -20 



-16 -14 



5 Discussion 



We have introduced general tools (the PRF and the ARF) to study the advance of both 
the phase and the amplitude variables for dynamical systems having a limit cycle 
attractor. These tools allow us to study variations of these variables under general 
perturbation hypotheses and extend the concept of infinitesimal PRCs which assumes 
the validity of the phase-reduction and is only true under strong hyperbolicity of the 
limit cycle or under weak perturbations. In fact, the PRFs and ARFs are first-order 
approximations of the actual variation of the phase and the amplitude, respectively, 
and so they are supposed to work mainly for weak perturbations; however, being an 
extension away from the limit cycle makes them more accurate than the PRCs even 
under strong perturbations. We thus claim that the phase-reduction has to be used 
with caution since assuming it by default may lead to completely wrong predictions 
in synchronization problems. We are not dismissing phase-reduction but trying to 
show the limits beyond which an extended scenario is required. 
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We have presented a computational analysis to understand the contribution of tran- 
sient effects in first-order predictions of the phase response, focusing on the impor- 
tance of the hyperbolicity of the limit cycle, but also on the relative positions of the 
isochrons with respect to the limit cycle. 

In the examples studied, subject to pulse-train stimuli, we have compared the pre- 
dictions obtained both with the new 2D map defined from the PRF and ARF and 
the ID map defined from the classical PRC. Using rotation numbers, we have shown 
differences up to two orders of magnitude in favor of the 2D predictions, especially 
when the stimulation frequency is high or the stimulus is too strong. These results 
confirm previous numerical experiments with specific oscillators; see [22]. On the 
other hand, we have found that both weak hyperbolicity of the limit cycle and "weak 
trans versality" of isochrons to it are important factors, although the role of hyperbol- 
icity seems to be more crucial. In this paper, these achievements have been tested in a 
canonical model allowing comparisons with the exact solutions, and other numerical 
tests have been applied in a conductance-based model. The technique can be ap- 
plied to other neuron models, and not necessarily for planar systems; n -dimensional 
systems would only require an additional computational difficulty in computing the 
associated (n — 1) ARFs. 

We would like to emphasize the importance of having good methods to compute 
isochrons (see [8-12]) since they are the cornerstone to study these transient phenom- 
ena that we have observed. They can be useful, not only for the problem illustrated 
here, but for other purposes like testing how far the experimentally recorded phase 
variations are from the theoretically predicted ones. In fact, they are the key con- 
cept to be able to predict the exact phase variation since, theoretically, if we know 
the parameterization K that gives the isochrons, the problem reduces to solving, at 
each step, (x, y) = K(0, a) and (V, y f ) = K(0 f , a'), where (x, y) is the point in the 
phase space where the pulse perturbation, ew, is applied and (x f , y f ) = (x,y) + ew. 
Indeed, the PRFs and ARFs can be computed knowing only the first order in K\ 
in principle, then they are valid only for weak perturbations, but easier to compute. 
Other refinements could be obtained by computing second order PRFs and ARFs by 
using the second-order approximations of the isochrons. Further extensions include 
also the possibility of computing response curves for long (in time) stimulus rather 
than pulsatile stimuli. 
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Appendix: The Vector Field for the A-Curves 

We prove here that given an analytic local diffeomorphism K, as in (5), satisfying 
(6), the A-curves are the orbits of a vector field Z, satisfying [X, Z] = [Z,X]= 0. 

This is equivalent to proving that DXZ = DZX. 

Taking derivatives with respect to 6 in Eq. (6), we get 

+ jvd^deK = (DX o K)doK, 

and using (14), we get 

By the chain rule, 

(DZ o K)(^d 0 + jad a ^K = (DX o K)(Z o K), 

and again, by the invariance equation (6), we obtain 

(DX o K)(Z oK) = (DZ o K)(X o K). (36) 
as we wanted to prove. 



References 

1. Izhikevich EM: Dynamical Systems in Neuroscience: The Geometry of Excitability and Bursting. 
Cambridge: MIT Press; 2007. [Computational Neuroscience.] 

2. Ermentrout B, Terman D: Mathematical Foundations of Neuroscience. New York: Springer; 2010. 

3. Schultheiss NW, Prinz AA, Butera RJ: Phase Response Curves in Neuroscience Theory, Experiment, 
and Analysis. Berlin: Springer; 2012. [Springer Series in Computational Neuroscience, vol 6.] 

4. Gutkin BS, Ermentrout GB, Reyes AD: Phase-response curves give the responses of neurons to 
transient inputs. J Neurophysiol 2005, 94(2):1623-1635. 

5. Ermentrout GB, Beverlin B, Troyer T, Netoff TI: The variance of phase-resetting curves. J Comput 
Neurosci 201 1, 31(2):185-197. 

6. Oh M, Matveev V: Non-weak inhibition and phase resetting at negative values of phase in 
cells with fast-slow dynamics at hyperpolarized potentials. J Comput Neurosci 2011, 31:31-42. 
doi: 1 0. 1 007/s 1 0827-0 1 0-0292-x. 

7. Winfree AT: Patterns of phase compromise in biological cycles. J Math Biol 1974/1975, 1:73-95. 

8. Guillamon A, Huguet G: A computational and geometric approach to phase resetting curves and 
surfaces. SIAMJAppl DynSyst 2009, 8(3): 1005-1042. doi: 10. 1137/080737666. 

9. Huguet G, de la Llave R: Computation of limit cycles and their isochrons: fast algorithms and 
their convergence. SIAM J Appl Dyn Syst 2013, in press. 

10. Osinga HM, Moehlis J: Continuation-based computation of global isochrons. SIAM J Appl Dyn 
Syst 2010, 9(4): 1201-1228. doi: 10. 1 137/090777244. [With online multimedia enhancements.] 



Springer 



Page 26 of 26 



O. Castejon et al. 



1 1 . Sherwood WE, Guckenheimer J: Dissecting the phase response of a model bursting neuron. SIAM 
JApplDyn Sys* 2010, 9(3):659-703. doi: 10.1 137/090773519. 

12. Mauroy A, Mezic I: On the use of Fourier averages to compute the global isochrons of 
(quasi)periodic dynamics. Chaos 2012, 22:033112. 

13. Ermentrout GB, Kopell N: Multiple pulse interactions and averaging in systems of coupled neural 
oscillators. J Math Biol 1991, 29(3):195-217. 

14. Brown E, Holmes P, Moehlis J: On the phase reduction and response dynamics of neural oscillator 
populations. Neural Comput 2004, 16:673-715. 

15. Oprisan SA, Canavier C: Stability analysis of rings of pulse-coupled oscillators: the effect of phase 
resetting in the second cycle after the pulse is important at synchrony and for long pulses. Differ 
Equ Dyn Syst 2001, 9:243-258. 

16. Oprisan SA, Prinz A A, Canavier CC: Phase resetting and phase locking in hybrid circuits of 
one model and one biological neuron. Biophys J 2004, 87(4):2283-2298. doi:10.1529/biophysj.l04. 
046193. 

17. Maran S, Canavier C: Using phase resetting to predict 1:1 and 2:2 locking in two neuron networks 
in which firing order is not always preserved. J Comput Neurosci 2008, 24:37-55. doi: 10. 1007/ 
sl0827-007-0040-z. 

18. Lin KK, Wedgwood KCA, Coombes S, Young LS: Limitations of perturbative techniques in the 
analysis of rhythms and oscillations. J Math Biol 2013, 66(1-2): 139-161 . doi:10.1007/s00285-012- 
0506-0. 

19. Ermentrout GB: Type I membranes, phase resetting curves, and synchrony. Neural Comput 1996, 
8:979-1001. 

20. Guckenheimer J: Isochrons and phaseless sets. J Math Biol 1974/1975, l(3):259-273. 

21. Cabre X, Fontich E, de la Llave R: The parameterization method for invariant manifolds. III. 
Overview and applications. J Differ Equ 2005, 218(2):444-515. 

22. Rabinovitch A, Friedman M: Fixed points of two-dimensional maps obtained under rapid stimu- 
lations. Phys Lett A 2006, 355(4-5):319-325. doi:10.1016/j.physleta.2006.02.059. 

23. Takeshita D, Feres R: Higher order approximation of isochrons. Nonlinearity 2010, 23(6): 1303- 
1323. doi:10.1088/0951-7715/23/6/004. 

24. Suvak O, Demir A: Quadratic approximations for the isochrons of oscillators: a general theory, 
advanced numerical methods and accurate phase computations. IEEE Trans Comput-Aided Des 
Integr Circuits Syst 2010, 29:1215-1228. 



4^ Springer 



